function CME = CMEstackFiltering(fileStations,OptsGen,OptStack,NC,Nstaz)
%CMEstackFiltering - Common mode error stack filtering
%
%   CME = CMEstackFiltering(fileStations,OptsGen,OptStack,Ncomp,Nstaz)
%
% This function carries out the CME filtering of the time series of the 
% GNSS stations listed in the Excel or ASCII file fileStations on the
% basis of the chosen stacking approach (see below about possible options).
% If data about detrended, zero mean time series are available in a tsData
% file, they can be used (see below about the input variable OptStack).
% It they are unavailable, or the user wants to recaculate them, the 
% detrended, zero-mean time series are computed and the corresponding 
% tsData files are upgraded.
% As the calculations are completed, the tsData files of involved stations 
% are automatically upgraded and the corresponding preprocessed mom files 
% are generated.
%
% The output CME is a struct variable whose fields are:
%       deltat : vector [t1 t2] of initial and final time (serial form)
%           Ivv: nt-by-ns matrix, where nt is the length of the time vector 
%                t1->t2 and ns is the number of stations, where Ivv(k,h) is 
%                true if at the time t(k) there are valid values for all
%                components of the station h
%          CMEE: nt-by-1 vector (options OptStack 1-3, 5-7, see below) or
%                nt-by-ns matrix (OptStack 4 or 8) of common mode error
%                for East component
%          CMEN: as above, North component
%          CMEV: as above, Vertical component 
% If the computations cannot be carried out, CME = [] is returned.
%
% Each station must be represented by its standard 4-character name in
% fileStation.  
% In the case of an Excel file, the station names must be placed in the 
% first column (no more than the first column is read). 
% In the case of an ASCII file, only a column is admitted. 
% If fileStations is undefined or empty, the filename can be interactively 
% managed.
%
% OptsGen is the struct variable with the general StaVel options, generated
% by means of geneOpts function. If OptsGen is undefined or empty, the name 
% of the file carrying such a struct variable can be interactively managed.
% In this case, the file must have a field whose name is either Opts or
% OptsGen.
%
% OptStack is the option on stacking. The admitted values are:
%   1:  simple stacking;
%   2:  stacking weighted by errors on daily positions;
%   3:  stacking weighted by distance of each station from the network 
%       center;
%   4:  stacking weighted by correlation between stations time series;
%   5:  simple stacking, with recalculation of possible existing 
%       detrended zero-mean time series;
%   6:  stacking weighted by errors on daily positions, with recalculation 
%       of possible existing detrended zero-mean time series;
%   7:  stacking weighted by distance of each station from the network 
%       center, with recalculation of possible existing detrended zero-mean
%       time series;
%   8:  stacking weighted by correlation between stations time series, with 
%       recalculation of possible existing detrended zero-mean time series;
%   9:  simple stacking with minimal detrending, i.e. detrending based on 
%       simple least square fitting to a straight line (in this case, the 
%       Hector-based trend is not used);
%   10.	stacking weighted by errors on daily positions with minimal 
%       detrending;
%   11.	stacking weighted by distance of each station from the network 
%       center with minimal detrending;
%   12.	stacking weighted by correlation between stations time series with 
%       minimal detrending.
% It OptStack is undefined, empty or is not an admitted value, it can be 
% interactively chosen.  
%
% Ncomp is the number of components. Allowed values: 2 (E-N), 3 (E-N-V).
% If Ncomp is undefined, empty or is not an allowed value, Ncomp = 3 is 
% used.
%
% Nstaz is the minimum number of stations to be used for stacking at each
% time. If, at a time t(k), less than Nstaz stations have data, CME(k) = 0 
% for each component. 
% Nstaz must be at least 3, i.e. at a time t(k), data from at least 3
% stations are necessary to allow the computation of CME(k). 
% If Nstaz is undefined, empty or lower than 3, the default Nstaz = 3 is 
% used.
% If Nstaz is Inf or higher than the number of stations, CMEs are non-zero 
% only at the times where all stations have data. 
%
% See also StaVelMain, detrendZeroMean, geneOpts, OptsGenLoad.

% G. Teza, 2022

if (nargin < 5) || isempty(Nstaz) || (numel(Nstaz) > 1) || (Nstaz < 3)
    Nstaz = 3;
end
if (nargin < 4) || isempty(NC) || (numel(NC) > 1) || ~ismember(NC,[2,3])
    NC = 3;
end
if NC == 2
    VC = [0 1];
else
    VC = [0 1 2];
end
if (nargin < 3) || (numel(OptStack) ~= 1) || ~ismember(OptStack,1:12)
    OptStack = [];
end
if nargin < 2
    OptsGen = [];
end
if nargin < 1
    fileStations = [];
end

b = readStationList(fileStations);
if isempty(b)
    CME = [];
    return
else
    ns = size(b,1);
    if ns < 3
        CME = [];
        return
    end
end

if isinf(Nstaz) || (Nstaz > ns)
    Nstaz = ns;
end

OptsGen = OptsGenLoad(OptsGen);

dirTs  = OptsGen.dirTs;
dirPreCME = OptsGen.dirPreCME;
AddTs = OptsGen.AddTs;
AddPreCME = OptsGen.AddPreCME;

if isempty(OptStack) 
    figIn = figure('Name','CME stack filtering');
    set(figIn,'Units','normalized','OuterPosition',[0 0 1 1]);
    imshow([])
    title('Common Mode Error filtering option','FontSize',24,'Color','b');  
    axis equal;
    lsgen = {...
        '1.  SIMPLE STACKING (AVAILABLE DETRENDED DATA)',...
        '2.  STACKING WEIGHTED BY ERRORS ON DAILY POSITIONS (AVAILABLE DETRENDED DATA)',...
        '3.  STACKING WEIGHTED BY DISTANCE OF EACH STATION FROM THE NETWORK CENTER (AVAILABLE DETRENDED DATA)',...
        '4.  STACKING WEIGHTED BY CORRELATION BETWEEN STATIONS TIME SERIES (AVAILABLE DETRENDED DATA)',...
        '5.  SIMPLE STACKING, WITH PRELIMINAR DETRENDING',...
        '6.  STACKING WEIGHTED BY ERRORS ON DAILY POSITIONS, WITH PRELIMINAR DETRENDING',...
        '7.  STACKING WEIGHTED BY DISTANCE OF EACH STATION FROM THE NETWORK CENTER, WITH PRELIMINAR DETRENDING',...
        '8.  STACKING WEIGHTED BY CORRELATION BETWEEN STATIONS TIME SERIES, WITH PRELIMINAR DETRENDING',...
        '9.  SIMPLE STACKING, MINIMAL DETRENDING',...
        '10. STACKING WEIGHTED BY ERRORS ON DAILY POSITIONS, MINIMAL DETRENDING',...
        '11. STACKING WEIGHTED BY DISTANCE OF EACH STATION FROM THE NETWORK CENTER, MINIMAL DETRENDING',...
        '12. STACKING WEIGHTED BY CORRELATION BETWEEN STATIONS TIME SERIES, MINIMAL DETRENDING'};
    lbIn = uicontrol(figIn,'Style','listbox',...
        'string',lsgen,...
        'FontSize',14,...
        'Units','normalized','Position', [0.1 0.3 0.85 0.40],...
        'Callback', 'uiresume(gcbf)');
    uiwait(figIn);
    OptStack = lbIn.Value;
    delete(lbIn);
    close(figIn);
end


% data read: 
t1 = zeros(1,ns);   % to define min t...
t2 = zeros(1,ns);   % ... and max t
Lat = nan(1,ns);
Lon = nan(1,ns);
MC = cell(1,ns);
nameS = repmat('    ',10,1);
Isok = false(1,ns);
for k = 1:ns
    
    namek = b{k,1};
    namek = strtrim(namek);     % to remove possible leading and trailing spaces
    namek = upper(namek);       % to warrant station name with capital letters
    if isempty(namek)
        continue
    end
    nameS(k,:) = namek;
        
    fileTsk = fullfile(dirTs,[namek AddTs '.mat']);
    tsk = tsDataFileIn(fileTsk);
    if isempty(tsk)
        continue
    end
    if isempty(tsk.detrendedZeroMeanE) || (OptStack >= 5)
        if OptStack >= 9
            [tsk,IT] = detrendZeroMeanMinimal(tsk,NC);
        else
            [tsk,IT] = detrendZeroMean(tsk,NC);
        end
        if IT
            ts = tsk; 
            save(fileTsk,'ts');
        else
            fprintf('Station %s - detrending impossible \n',namek);
            continue
        end
    end
    Isok(k) = true;
    tink = tsk.t;
    nink = numel(tink);
    sEink = tsk.sE;
    sNink = tsk.sN;
    sVink = tsk.sV;
    Lat(k) = tsk.statLat;
    Lon(k) = tsk.statLon;
    %
    tEk = tsk.detrendedZeroMeanE.t;
    oEk = tsk.estimatedTrendE.o;        % observation
    dEk = tsk.detrendedZeroMeanE.d;     % detrended
    IEk = ismember(tink,tEk);
    %
    tNk = tsk.detrendedZeroMeanN.t;
    oNk = tsk.estimatedTrendN.o;
    dNk = tsk.detrendedZeroMeanN.d;
    INk = ismember(tink,tNk);
    % 
    if NC == 3
        tVk = tsk.detrendedZeroMeanV.t;
        oVk = tsk.estimatedTrendV.o;
        dVk = tsk.detrendedZeroMeanV.d;
        IVk = ismember(tink,tVk);
    else
        IVk = true(nink,1);
    end
    Ik = logical(IEk.*INk.*IVk);
    if (NC == 2) || isempty(tsk.detrendedZeroMeanV)
        Mk = zeros(nink,7);
    else
        Mk = zeros(nink,10);
    end
    Mk(:,1) = tink;
    Mk(IEk,2) = oEk;
    Mk(IEk,3) = dEk;
    Mk(:,4) = sEink;
    Mk(INk,5) = oNk;
    Mk(INk,6) = dNk;
    Mk(:,7) = sNink;
    if (NC == 3) && ~isempty(tsk.detrendedZeroMeanV)
        Mk(IVk,8) = oVk;
        Mk(IVk,9) = dVk;
        Mk(:,10) = sVink;
    end
    Mk = Mk(Ik,:);
    MC(k) = {Mk};
    t1(k) = min(Mk(:,1));
    t2(k) = max(Mk(:,1));
end

% Generation of correlated matrices, taking into account problems of 
% possible data missing:
t1t = min(t1);
t2t = max(t2);
t = (t1t:t2t)';
nt = numel(t);

oE = nan(nt,ns);
dE = nan(nt,ns);
sE = nan(nt,ns);
oN = nan(nt,ns);
dN = nan(nt,ns);
sN = nan(nt,ns);
Ivv = false(nt,ns); % indices of valid values. Ivv(k,h) is 1 if for t(k)
                    % all the components for the station h are valid 
if NC == 3
    oV = nan(nt,ns);
    dV = nan(nt,ns);
    sV = nan(nt,ns);
end
for k = 1:ns
    if Isok(k)
        Mk = MC{k};
        tk = Mk(:,1);
        Ivvk = ismember(t,tk);
        Ivv(:,k) = Ivvk;
        oE(Ivvk,k) = Mk(:,2);
        dE(Ivvk,k) = Mk(:,3);
        sE(Ivvk,k) = Mk(:,4);
        oN(Ivvk,k) = Mk(:,5);
        dN(Ivvk,k) = Mk(:,6);
        sN(Ivvk,k) = Mk(:,7);
        if NC == 3
            oV(Ivvk,k) = Mk(:,8);
            dV(Ivvk,k) = Mk(:,9);
            sV(Ivvk,k) = Mk(:,10);
        end
    end
end

if OptStack >= 9
    OptStack = OptStack-8;
elseif OptStack >= 5
    OptStack = OptStack-4;
end

if ismember(OptStack,1:3)
    if OptStack == 1
        % Simple stacking
        CMEE = nansum(dE,2)./nansum(Ivv,2); %#ok<*NANSUM>
        CMEN = nansum(dN,2)./nansum(Ivv,2);
        if NC == 3
            CMEV = nansum(dV,2)./nansum(Ivv,2);
        end
        strStack = 'Stacking';
    elseif OptStack == 2
        % Stacking weighted with standard deviation
        wE = 1./(sE.^2);
        CMEE = nansum(dE.*wE,2)./nansum(wE,2);
        wN = 1./(sN.^2);
        CMEN = nansum(dN.*wN,2)./nansum(wN,2);
        if NC == 3
            wV = 1./(sV.^2);
            CMEV = nansum(dE.*wV,2)./nansum(wV,2);
        end
        strStack = 'Standard deviation weighted stacking';
    else
        % Stacking weighted with distance
        x = nan(1,ns);
        y = nan(1,ns);
        kmin = find(Isok,1);    % first valid station
        for k = 1:ns
            if Isok(k)
                Latk = Lat(k);
                Lonk = Lon(k);
                if k == kmin
                    [xk,yk,utmzone1,utmhemi1] = wgs2utm(Latk,Lonk);
                else
                    [xk,yk,~,~] = wgs2utm(Latk,Lonk,utmzone1,utmhemi1);
                end
                x(k) = xk;
                y(k) = yk;
            end
        end
        x0 = nanmean(x); %#ok<*NANMEAN>
        y0 = nanmean(y);
        sqDistCM = (x-x0).^2+(y-y0).^2;
        w = 1./sqDistCM;
        CMEE = nansum(dE.*w,2)./nansum(w,2);
        CMEN = nansum(dN.*w,2)./nansum(w,2);
        if NC == 3
            CMEV = nansum(dV.*w,2)./nansum(w,2);
        end
        strStack = 'Distance-weighted stacking';
    end
elseif OptStack == 4
    % initialization:
    CMEE = zeros(nt,ns);
    CMEN = zeros(nt,ns);
    rE = zeros(ns);
    rN = zeros(ns);
    wE = 1./(sE.^2);
    wN = 1./(sN.^2);
    if NC == 3
        CMEV = zeros(nt,ns);
        rV = zeros(ns);
        wV = 1./(sV.^2);
    end
    % correlation between stations:
    for l = 1:ns
        for m = 1:ns
            rElm = corr([dE(:,l),dE(:,m)],'rows','complete');
            rE(l,m) = rElm(1,2);
            rNlm = corr([dN(:,l),dN(:,m)],'rows','complete');
            rN(l,m) = rNlm(1,2);
            if NC == 3
                rVlm = corr([dV(:,l),dV(:,m)],'rows','complete');
                rV(l,m) = rVlm(1,2);
            end
        end
    end
    % computation of CMEs:
    for k = 1:ns
        % E:
        numEk = nan(nt,ns);
        denEk = nan(nt,ns);
        for h = 1:ns
            numEkh = dE(:,h).*wE(:,h)*rE(k,h);
            denEkh = wE(:,h)*rE(k,h);
            numEk(:,h) = numEkh;
            denEk(:,h) = denEkh;
        end
        CMEEk = nansum(numEk,2)./nansum(denEk,2);
        CMEE(:,k) = CMEEk;
        % N:
        numNk = nan(nt,ns);
        denNk = nan(nt,ns);
        for h = 1:ns
            numNkh = dN(:,h).*wN(:,h)*rN(k,h);
            denNkh = wN(:,h)*rN(k,h);
            numNk(:,h) = numNkh;
            denNk(:,h) = denNkh;
        end
        CMENk = nansum(numNk,2)./nansum(denNk,2);
        CMEN(:,k) = CMENk;
        % V:
        if NC == 3
            numVk = nan(nt,ns);
            denVk = nan(nt,ns);
            for h = 1:ns
                numVkh = dV(:,h).*wV(:,h)*rV(k,h);
                denVkh = wV(:,h)*rV(k,h);
                numVk(:,h) = numVkh;
                denVk(:,h) = denVkh;
            end
            CMEVk = nansum(numVk,2)./nansum(denVk,2);
            CMEV(:,k) = CMEVk;
        end
    end
    strStack = 'Correlation-based stacking';
end

% consistency on minimum number of good stations:
Inov = sum(Ivv,2) < Nstaz;  % non valid times for CMEs 
CMEE(Inov,:) = 0;
CMEN(Inov,:) = 0;
if NC == 3
    CMEV(Inov,:) = 0;
end
% suppression of points having too high CME:
if NC == 2
    IexCME = (abs(CMEE)>3*mean(abs(CMEE)))|(abs(CMEN)>3*mean(abs(CMEN)));
else
    IexCME = (abs(CMEE)>3*mean(abs(CMEE)))|(abs(CMEN)>3*mean(abs(CMEN)))|(abs(CMEV)>3*mean(abs(CMEV)));
    CMEV(IexCME) = 0;
end
CMEE(IexCME) = 0;
CMEN(IexCME) = 0;

% filtering:
for k = 1:ns
    if Isok(k)
        % simple Filtering:
        oEF = oE-CMEE;   % (CMEs are column vectors or ns-column matrices)
        dEF = dE-CMEE;
        oNF = oN-CMEN;
        dNF = dN-CMEN;
        if NC == 3
            oVF = oV-CMEV;
            dVF = dV-CMEV;
        end
        % reconstruction of original time series
        Ivvk = Ivv(:,k);
        tk = t(Ivvk);
        oEk = oE(Ivvk,k);
        dEk = dE(Ivvk,k);
        oEFk = oEF(Ivvk,k);
        dEFk = dEF(Ivvk,k);
        oNk = oN(Ivvk,k);
        dNk = dN(Ivvk,k);
        oNFk = oNF(Ivvk,k);
        dNFk = dNF(Ivvk,k);
        if NC == 3
            oVk = oV(Ivvk,k);
            dVk = dV(Ivvk,k);
            oVFk = oVF(Ivvk,k);
            dVFk = dVF(Ivvk,k);
        end
        % For the generation of output  files and tsData file upload
        namek = nameS(k,:);
        fileTsk = fullfile(dirTs,[namek AddTs '.mat']);
        tsk = tsDataFileIn(fileTsk);
        tsk.CMEfiltered.t = tk;
        tsk.CMEfiltered.E = oEFk;
        tsk.CMEfiltered.N = oNFk;
        if NC == 3
            tsk.CMEfiltered.V = oVFk;
        else
            tsk.CMEfiltered.V = [];
        end
        tsk.CMEfiltered.Method = strStack;
        ts = tsk;
        save(fileTsk,'ts');
        MosMomk = MosMomGen(tsk,1);
        Mosk = MosMom2MosNeu(MosMomk);
        nameOutMomk = fullfile(dirPreCME,[namek AddPreCME]);
        writeNeuMom(nameOutMomk,VC,tsk,Mosk,'PreCME');
        % figures:
        figure('name',sprintf('%s STATION - %s',namek,strStack),...
            'units','normalized','outerposition',[0 0 1 1]);
        %
        subplot(NC,2,1)
        plot(tk,oEk,'.k');
        hold on; 
        plot(tk,oEFk,'.r');
        legend('original ts','CM-filtered');
        datetick('x',28); grid on; 
        ylabel('\itEast \rm(m)');
        %
        subplot(NC,2,2)
        plot(tk,dEk,'.b');
        hold on; 
        plot(tk,dEFk,'.g');
        legend('original detrended zero-mean ts','CM-filtered detrended zero-mean ts');
        datetick('x',28); grid on; 
        ylabel('\itEast \rm(m)');
        %
        subplot(NC,2,3)
        plot(tk,oNk,'.k');
        hold on; 
        plot(tk,oNFk,'.r');
        legend('original ts','CM-filtered');
        datetick('x',28); grid on; 
        ylabel('\itNorth \rm(m)');
        %
        subplot(NC,2,4)
        plot(tk,dNk,'.b');
        hold on; 
        plot(tk,dNFk,'.g');
        legend('original detrended zero-mean ts','CM-filtered detrended zero-mean ts');
        datetick('x',28); grid on; 
        ylabel('\itNorth \rm(m)');
        %
        if NC == 3
            subplot(NC,2,5)
            plot(tk,oVk,'.k');
            hold on; 
            plot(tk,oVFk,'.r');
            legend('original ts','CM-filtered');
            datetick('x',28); grid on; 
            ylabel('\itVertical \rm(m)');
            %
            subplot(NC,2,6)
            plot(tk,dVk,'.b');
            hold on; 
            plot(tk,dVFk,'.g');
            legend('original detrended zero-mean ts','CM-filtered detrended zero-mean ts');
            datetick('x',28); grid on; 
            ylabel('\itVertical \rm(m)');
        end
    end
end

CME = struct('CMEE',CMEE);
CME.CMEN = CMEN;
if NC == 3
    CME.CMEV = CMEV;
else
    CME.CMEV = [];
end
CME.deltat = [t(1) t(end)];
CME.Ivv = Ivv;